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1  Abstract 


In  this  project  we  aim  to  construct  a  high  fidelity  boundary  condition  module  for  Maxwell’s  equations 
that  can  be  interfaced  with  major  time-domain  electromagnetics  solver  systems.  There  is  ample  need 
in  the  EM  modeling  community  for  reliable  and  stable  far  field  boundary  conditions  of  high  accuracy. 
Most  existing  methods  are  limited  in  one  or  more  of  these  requirements,  and  recent  developments  in 
the  CRBC  procedure  (as  originally  presented  by  Hagstrom  and  Warburton  in  2009),  have  made  the 
technique  an  attractive  candidate  for  implementation  in  multi-purpose  solvers.  In  phase-I  of  this  project 
we  implemented  and  improved  upon  many  aspects  of  this  method,  particularly  in  light  of  the  needs  of 
high  order  accurate  Maxwell  equations  solvers  (based  on  the  discontinuous  Galerkin  method).  Error 
bounds  were  computed  and  demonstrated  for  a  number  of  cases.  We  continue  in  the  second  phase  of 
this  project  to  improve  upon  the  robustness  of  this  method,  as  we  develop  a  software  platform  which 
shall  be  its  flagship  (and  open  source)  implementation.  In  this  first  quarterly  report  we  present  initial 
progress  to  this  end,  showing  recent  mathematical  developments,  project  coordination  plans  and  some 
results  from  our  initial  experience  with  MEEP. 
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2  Introduction 


In  this  project,  HyPerComp  is  collaborating  with  Prof.  Thomas  Hagstrom  and  his  research  group  at 
the  Southern  Methodist  University  (SMU).  Roles  of  the  two  organizations  are  very  broadly  divided 
into  mathematical  method  development  (led  by  SMU)  and  implementation,  software  development  and 
maturation  (led  by  HyPerComp).  The  project  is  coordinated  via  a  series  of  in-person  and  telephone 
meetings.  As  of  May  17,  2013  we  have  been  conducting  weekly  telephone  meetings.  Two  students,  John 
Lagrone  and  Fritz  Juhnke  have  been  included  in  the  team  and  have  been  actively  participating  in  the 
work  so  far. 

Tasks:  The  following  is  a  list  of  tasks  to  be  performed  in  this  project. 

1.  Project  Formulation 

2.  Software  Development 

3.  Verification  Validation 

4.  Coupling 

5.  Efficiency  Testing 

6.  Release  of  software 

7.  Documentation 

8.  Sustainability  Plan 

9.  User  Support 

At  present,  we  are  working  on  a  refined  formulation  of  the  CRBC  method  and  testing  it  in  sample 
problems.  Primary  concerns  pertaining  to  method  stability  at  corners,  particularly  in  3D  are  being 
addressed.  CRBC  implementations  in  finite  difference  schemes,  DC  (in  FORTRAN  as  well  as  in  MAT- 
LAB)  are  available  from  prior  research  in  this  project,  for  testing.  We  hope  to  host  a  telephone  meeting 
with  our  project  manager  from  the  Army  in  July,  coinciding  with  Prof.  Hagstrom’s  visit  to  HyPerComp. 

We  are  presently  aiming  to  integrate  the  CRBC  module  with  the  following  codes: 

•  HDphysics  from  HyPerComp,  a  high  order  DC  based  solver 

•  MEEP  from  MIT,  an  open  source  FDTD  code 

•  cgmx  part  of  “Overture”  suite  of  simulation  codes  from  LLNL  -  high  order  finite  differences, 
second  order  PDEs 

•  CLAWPACK  a  finite  difference  suite  of  solvers  from  U. Washington 

Students  from  SMU  shall  initially  focus  on  an  FDTD  implementation  with  MEEP  and  begin  to 
identify  software  needs.  We  are  in  the  process  of  developing  software  requirements  for  each  of  the  systems 
mentioned  above,  so  that  we  can  outline  a  common  implementation  of  the  method  and  programming 
techniques.  This  shall  be  discussed  in  the  forthcoming  report. 
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3  CRBCs  and  Corner  Conditions  for  the  TM  Maxwell  System 


Consider  the  TM  Maxwell  system: 


dH^  1  dE^ 

dt  fj,  dy 

any  i  dE^ 

dt  fjL  dx 

dE^  1  any  i  an^ 

at  e  dx  e  dy 


(1) 

(2) 

(3) 


and  set  c  =  —)= . 


CRBC  on  an  arbitrary  edge 


Consider  a  portion  of  the  radiation  boundary  with  unit  normal  n  pointing  outward  from  the  compu¬ 
tational  domain  and  unit  tangent  vector  t: 


n  = 


(4) 


(Notice  we  have  chosen  a  specific  orientation  for  simplicity.)  Introducing  angles  cj)j, 

d 


j  =  l,...,P,we 


rewrite  the  Maxwell  system  in  terms  of  normal  and  tangential  derivatives,  ^  and  ^  and  recursively 
replace  normal  derivatives  with  interpolation  operators.  For  outgoing  waves 


a  cos  (pj  a  sin^  pj 
dn  c  at  cT  cos  pj  ’ 


(5) 


and  for  incoming  waves 


a  cos  Pi  a  sin^  Pi  ,  . 

dn  c  at  cT  cos  pj 

Although  not  necessary,  it  is  convenient  to  recast  the  first  two  equations  using  the  normal  and  tangential 
components  of  H: 


+  nyHy)  + 

+n,Hy)- 


1  dE^ 
y  dr 

1  as^ 

y  dn 


0 

0 


(7) 

(8) 


Introducing  auxiliary  variables  F/J,  iJj,  and  Ej  with  Hq,  and  Eg  coinciding  with  the  trace  of  the 
solution  at  the  boundary  (or  in  the  DC  context  the  boundary  state)  we  solve  for  j  =  1, . . .  P 


¥t 


I  TJX  rrV  \  COS  Pi  aE^_^ 


sin^  pj 


yc  dt  ficT  cos  pj 


E: 


i-i 


TTX  ttv\  cos  pjdE^j 


1  sin  Pj^^ 
ycT  cos  pj  ^ 


(9) 
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and  for  j  =  0, 


dE^-i  cos4>j  d 


dt 


1  sin  (j)^ 
ecT  cos  (pj 


dt 
1  d 


(  nyHj_i+nxHyi) 


{-UyHl,  +  nxHy)  +  -  —  {nxHl,  +  riyRy)  = 


dr 

dE^  cos  4>j  d  ,  y. 

^  {-nyH^  +  rixHy) 


dt 


ec 


dt 


■  ecT  cost  ^  ^  +  ^y^j)  ’ 


,P 


-(uxH^^nyHy) 


dt 


1  dE^ 
^  dr 


=  0. 


(10) 


(11) 


These  are  3P  +  1  equations  for  3P  +  3  unknowns.  Two  additional  equations  are  obtained  first  by  using 
data  from  the  interior  for  outgoing  characteristics  and  terminating  the  incoming  characteristic  recursion. 
The  system  simplifies  if  we  introduce  the  normal  characteristic  variables: 


R^,y  =  E^J  ±  (-UyHj  +  UxB])  ,  =  UxH^  +  nyH^ 


then  the  recursions  take  the  form 


(1  +  COSpj) 


dR+  1  sin^ 


dt 


+  — - 


T  cospj 


R+,j-i  +  “  ■ 


IdH 


n,i-l 


(1  -  COSpj) 


-  ,  dR+j  1  sin^  (pj 


dt  T  cos  (pj 


R+j  +  -  - 


dr 

IdHr,. 


dr 


(12) 


(13) 


We  then  take 


dR^j-i  Isin^  (pj  IdHnj.i 

(l-cos((.j-)  + -- 


dt 


T  cos  (pj 


e  dr 


-  dR^j  Is\T?(pj  ldHn,j 

( 1  +  COS  (bj )  — — - H - = —  R-  7  + - — — . 

^  dt  T  cos(pj  e  dr 


dR^ 


ldR^,o\ 

\  dt  ) 


dt 


interior 


dR 

and  solve  (14)  in  increasing  j  for  .  The  termination  condition  is 

R+,p  =  0, 

which  allows  us  to  solve  (13)  in  decreasing  j  for  ■ 


(14) 


(15) 


(16) 


Corner 

We  now  consider  the  case  of  two  artificial  boundaries  meeting  at  a  corner  point.  Let  =  (u^juJ)^  be 
the  unit  outward  normal  for  the  first  boundary  and  =  (^^,71^)^  be  the  unit  outward  normal  for  the 
second.  Introduce 

7  =  7^  <  1. 
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Notice  that  if  the  angle  between  the  edges  is  obtuse  then  7  >  0.  Also  introduce  Sj  =  ±1,  j  =  1,2  with 
Sj  =  1  if  the  edge  orientation  is  into  the  corner  and  Sj  =  -1  otherwise.  We  now  solve  for  a  doubly- indexed 
array  of  auxiliary  variables  Hjf,  and  which  coincide  with  the  corner  values  of  the  auxiliary 

variables  E^’^ ,  7  when  j  =  0  and  Ej’^,  when  k  =  0.  Thus  the  auxiliary  variables 

with  index  0, 0  correspond  to  the  corner  values  of  the  actual  fields.  A  complete  set  of  equations  for  the 
corner  variables  is  obtained  by  applying  the  interpolation  conditions  to  replace  all  space  derivatives. 
We  do  so  by  solving  for  \JW  in  terms  of  n}  ■  YW  and  n?  ■  YW  where  W  is  any  function.  This  yields 

VIT  =  -VW)  (17) 

+  (1  -  -  7n^)(n^  •  YW) 


To  further  simplify  the  formula  we  use  the  fact  that  'v}'  =  'n^'  =  1  and  find  that 


=  1  -  (n;J,n^  + 

=  +  {nlf)  {{nlf  +  {nlf)  - 

=  (nlf  (nlf  +  {nlf  {nlf  -  2nlnlnlnl 
=  S'2 


where  we  have  introduced 


Then 


Similarly 


Thus 


o  12  12 

S  —  n^Uy  HyTL^. 


{n}-^v?)  = 


nl-nl{nlf  -nlnln 

n 


y"^y"^x  \ 

1  _  I 


Vy  '^X  ^X' 

nl{nlf  -  n'ln'ln: 

nlnlnl 


n:,{nl) 

,2 


=  5 


-y\ 

( ) 


n:,{nly 

1 


(l-7^)-i(n^-7ni)  =  - 


( 1 ) 


JL  ^  f(  2j9__  1  d  \ 

dx  S  drC  dri^  ) 

d  1  /  2  ^  1  ^  ) 


(18) 

(19) 


(20) 

(21) 

(22) 

(23) 


Then  all  derivatives  can  be  replaced  by  the  interpolation  conditions  (5)-(6)  associated  with  one  of 
the  edges.  We  also  assume  for  simplicity  that  we  are  using  the  same  angle  parameters  and  orders  on 
each  edge.  This  is  not  necessary,  but  it  is  what  we  usually  do  and  assuming  it  simplifies  the  notation. 


The  recursions  corresponding  to  (9)  are  the  easiest  to  write  down.  We  solve  for  j  = 
k  =  0,...,P 


cy 

Ft 


COS 


fiC 


dt 


A/'_„2rra:  .  ttV  \  COS(j>jdE-j, 


1  siF  (l)j 
ficT  cos(j)j 
1  sin^^j- 
ficT  cos(j)j 


(24) 
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and  for  j  =  0, . . . , P,  k  =  1, . . .  ,P 


d 


^  +  ntH 

ot 


1  ttV 


_9 

dt 


COS 


4>k^^lk-i  1  sin^cpk 


lie 


dt  iicT  cos  (f>k 


E- 


cosepk 


1  sin  dk 


lie  dt  iieT  cosdk 


Elk- 


(25) 


Lastly  from  the  (10)  we  derive  a  relationship  involving  both  recursions  which  we  impose  for  j,k  = 
1, . . . ,  P.  Using  (22)-(23)  to  express  the  space  derivatives  using  the  normals  we  hnd 


^dE^  odRy  .dRy  ^dH^  .dH^  „ 
eS^^-n;^—  +  M^—-ni^^-kni^^  =  0. 
"  dM  "  dn^ 


dt  y  dM  y  drp  ^  drd 
Now  replacing  the  normal  derivatives  by  (5)- (6)  we  derive: 


dn? 


(26) 


ec  dt 

sin^  (h^ 

+ - ^ — 

ecT  cos  (l>j 


+  El^_^ 


cos  4>j  d  I  I  ( 


■-Elk-i) 


+ 


'  Tjy 


1 - ^— 

ecT  cos  (f>j 

,  COS(j)k  d  /  2  (  jxx  TJX  t  ,  ^2  (  r 

*ecTcos<Pt  -  Kt-i)  *  "J  (U-i.i-1 ' 

u-HUd-<(CM- 


-EUd 

)) 
)) 


Ry,  , 


Ry,  , 


+ 


ecT cos  4>k  ~ 


)) 


"Ud) 


=  0 


(27) 


We  have  now  written  down  2P(P+ 1)  +  P^  =  3P^  +  2P  equations  in  3(P+ 1)^  =  3P^  +  6P  +  3  variables. 
Thus  4P  +  3  additional  equations  are  required.  We  first  incorporate  incoming  data  from  the  edges.  Here 
we  define  outgoing  characteristics  (ingoing  to  the  corner)  in  the  tangential  directions  for  each  edge. 
Then  for  fc  =  0, . . . ,  P  set 


Sk  -  Eo^k  +  sn 


^0,k 


•A 


(28) 


Then  with  the  correspondence 


H. 


0,k’> 


Elk) 


(  JTZ  ttX 


.Hi) 


edgej 


(29) 


we  can  impose 


dt  \  j 


(30) 
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Similarly  for  j  =  0, . . . ,  P  set 

^3  -  ^lo  +  S2- 

^/{<Hl,PnlHl,), 

(31) 

and  note  the  correspondence 

{El,,  HI,,  Hi 

o)-{El,Hl 

^j/^edge. 

(32) 

Then  we  can  impose 

dSl 

2 X  edge. 

3 

dt 

(33) 

Finally  we  can  impose  the  termination  conditions 

^+,3,P 

=  0  i  =  o,... 

,P 

(34) 

R+.P.k 

o' 

II 

o 

II 

.,P.. 

(35) 

These  constitute  4(P  +  1)  =  4P  +  4  equations,  one  more  than  we  can  use.  We  remove  an  equation  by 
only  imposing  the  sum  of  (30)  and  (33)  for  j  =  k  =  0. 


4  Meep 


iVfeep-4..2  and  related  supporting  software,  libctl-3.2. 1,  harminv,  harminv-develhave  been  installed 
in  Linux  systems.  Single  CPU  build  has  been  tested  and  verified  against  the  examples  posted  on  the 
Meep  website.  Parallel  version  has  not  been  tested  yet. 

In  Fedora  OS  systems,  default  configure  for  libctl-3.2.1  will  not  work,  instead,  we  need  to  use 
./configure  LDFLAGS—-lm 
to  configure  the  software. 

Packages  harminv  and  harminv-devel  will  be  needed  for  solutions’  Fourier  transforms  and  need 
to  be  installed  before  compile  Meep. 

Default  configure  for  Meep- 1.2  also  failed  in  our  Linux  system.  It  cannot  automatically  include  the 
Lapack  library  when  trying  to  link  the  software.  We  need  to  explicitly  add  the  Lapack  library  as  follows: 
./configure  LDFLAGS—-llapack  prefix— $Meep_home 
with  $  Meep  home  the  destination  of  Meep  installation. 

Post-processing  utility  hSutils- 1.12.1  is  also  installed.  There  are  a  lot  of  utilities  included  in  this 
package.  To  manipulate  unsteady  solution  easier,  we  developed  a  python  script  to  call  hStovtk  utility 
and  create  a  series  of  .vtk  file  from  the  big  .h5  file.  The  script  is  list  below 

# ! /usr/bin/python 
import  commands 
import  string 
import  sys 

argc  =  lenCsys . argv) 
if  argc  ! =  2 : 

print  "usage  error,  wrong  number  of  arguments" 

print  "usage  example:  . /convHSToVtk.py  sq_scatter-hz.h.5" 
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exit (-1) 

#get  the  final  time  step 
cmd  =  "h51s  "  +  sys . argv [argc-1] 
hBlso  =  commands .getoutput (cmd) 
hBlso  =  hBlso. replace("/Inf}" ,  "") 

Ispace  =  hBlso. rfindC  ") 

tend  =  int (hBlso [lspace+1 :] ) 

print  "final  time  step  is:  ",  tend-1 

fileRoot  =  sys . argv [argc-1] 
fileRoot  =  fileRoot [: len(fileRoot)-3] 
defaultVtk  =  f ileRoott" . vtk" 
print  fileRoot 

for  i  in  range(tend): 
bigNum  =  1000000+i 

newVtk  =  fileRoot+"_"+(str (bigNum) ) [1 :]+". vtk" 

cmdl  =  "hBtovtk  -t  "  +  str(i)  +  "  "  +  sys . argv [argc-1] 

commands . getoutput (cmdl) 

cmd2  =  "mv  "  +  defaultVtk  +  "  "  +  newVtk 

commands . getoutput (cmd2) 

print  newVtk  +  "  is  ready" 

With  .vtk  files,  we  can  easily  use  paraview  to  do  more  detailed  analysis. 

Meep  code  was  verified  with  several  examples  using  the  input  data  provided  in  the  Meep  Tutorial 
web  page.  Grid  resolution  and  pml-layer  thickness  have  been  studied  to  check  the  code  accuracy  and 
pml-layer  sensitivity. 

Figure  1  gives  the  solution  of  E^  2-D  held  in  a  a  straight  waveguide.  The  waveguide  has  e  =  12. 
Figure  2  depicts  the  solution  of  a  90°  bend  waveguide  solution  at  t  =  79.8, 139.8  and  199.8  respectively. 
To  further  check  the  code  capability,  a  3D  case  with  a  perfect-electric-conductor  sphere  put  at  the  center 
of  the  domain  is  calculated.  The  wavelength  is  set  to  A  =  Im,  and  the  source  frequency  nondimensional 
frequency  /  =  1  (which  indicates  the  real  frequency  =  f*c/X  =  3x  10®  Hz).  The  domain  is  set  to 
10  X  10  X  10,  and  the  radius  of  the  sphere  is  set  to  r  =  10A/(27r)  =  1.59m.  Figure  3  shows  the  E^  held 
distribution  at  non-dimensional  t  =  200.  For  this  run,  with  resolution=20,  it  takes  about  2  hours  in  an 
Intel  core-i7  CPU. 

In  Figure  2  ,  we  computed  the  held  patterns  for  light  propagating  around  a  waveguide  bend.  The 
results  are  only  visually  pretty  and  not  quantitatively  satisfying.  We’d  like  to  know  exactly  how  much 
power  makes  it  around  the  bend,  how  much  is  rehected,  and  how  much  is  radiated  away.  This  is  done  by 
running  the  example  case  bend-flux  in  the  Meep  example  website.  The  domain  size  is  16  x  32,  including 
the  pml-layers,  which  we  set  the  thickness  to  be  1.  We  follow  the  instruction  in  the  tutorial,  and  run 
with  resolution  10,  20  and  40.  Figure  4  depicts  the  comparison  of  the  case  with  different  resolution.  It 
can  be  seen  a  minimum  resolution  of  20  is  needed. 

In  addition,  we  tried  different  pml  thickness  with  resolution  20.  Figure  5  depicts  the  comparison  of 
the  case  with  different  pml  thickness.  This  is  to  check  the  robustness  of  the  radiation  BC  in  Meep  code 
and  see  how  much  we  can  improve  by  using  CRBC  module  in  this  project  in  the  future.  Obviously,  with 
very  thin  pml,  the  results  oscillate,  which  is  in  our  expectation.  However,  when  pml  thickness  increased 
to  2,  the  solution  is  also  far  away  from  the  thickness  1  benchmark,  which  we  may  need  to  find  out  later. 
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Figure  1:  Electric-field  component  E^  distribution  of  straight  waveguide  field  at  t  =  200 


Figure  2:  Electric-field  component  E^  distribution  at  of  the  90°  waveguide  field  at  t  =  79.8, 139.8  and 
199.8. 
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Figure  3:  Electric-field  component  E^  distribution  of  a  3D  case  at  t  =  200 
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Transmission  spectrum  around  a  waveguide  bend 


Transmission  spectrum  around  a  waveguide  bend 


Figure  4:  Transmission  spectrum  around  a  waveguide  bend  using  different  grid  resolution. 
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Transmission  spectrum  around  a  waveguide  bend 


Transmission  spectrum  around  a  waveguide  bend 


Figure  5:  Transmission  spectrum  around  a  waveguide  bend  using  different  pml  thickness. 
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